Maximally-localized Wannier functions for disordered systems: application to 

amorphous silicon 
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We use the maximally-localized Wannier function method to study bonding properties in amorphous 
silicon. This study represents, to our knowledge, the first application of the Wannier-function 
analysis to a disordered system. Our results show that, in the presence of disorder, this method 
is extremely helpful in providing an unambiguous picture of the bond distribution. In particular, 
defect configurations can be studied and characterized with a novel degree of accuracy that was not 
available before. 
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I. INTRODUCTION 

Since their introduction in 1937 Wannier functions!!] 
have played an important role in the theoretical study of 
the properties of periodic solids. Moreover the represen- 
tation of the electronic ground state of periodic systems 
in terms of localized Wannier or Wannier-like orbitals 
has recently attracted considerable attention due to the 
development of "order-N" methodsci and to the formu- 
lation of the modern theory of electronic polarizationB. 
In the case of finite systems, localized orbitals are widely 
used to describe and understand chemical concepts such 
as bonds, lone-pair orbitals, and .valence-electron charge 
distributions; different criteriaQcl have been developed 
for producing optimum localized orbitals. 

In periodic systems the determination of localized or- 
bitals or Wannier functions is,-much less trivialQ'El. Re- 
cently, Marzari and VanderbiliH developed a very practi- 
cal method for generating maximally-localized Wannier 
functions starting from the knowledge of the occupied 
Bloch states. This amounts to the generalization fox. 
periodic systems of the Boys' localized-orbital methods 
that is commonly used in quantum chemistry. The new 
technique has been successfully applied to crystal sys- 
tems and small moleculesu. Here we apply for the first 
time the same procedure to a disordered system, namely 
amorphous silicon. We show that, also in this case, the 
Wannier functions are extremely useful in providing a 
clear description of the relevant electronic and bonding 
properties. They also help in eliminating many of the 
ambiguities that are usually associated with identifying 
defects in a disordered system. 



II. METHOD 

The Wannier functions!*! are defined in terms of a uni- 
tary transformation of the occupied Bloch orbitals. Even 
for the case of a single band, the Wannier functions are 
not uniquely defined, due to the arbitrary freedom in the 



phases of the Bloch orbitals. In the multiband case this 
freedom becomes more general, and includes the choice 
of arbitrary unitary transformations among all the occu- 
pied orbitals at every paint in the Brillouin zone (BZ). 
Marzari and Vanderbiltfl resolve this indeterminacy by 
requiring that the total spread of the Wannier functions 
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be minimized in real space, in analogy with Boys' 
criteriontj for finite systems. In Eq. (jl) (...)„ indicates 
the expectation value with respect to t jjs n-th Wannier 
function w n (r). Marzari and Vanderbilttl have discussed 
how to properly define the operators r and r 2 in a peri- 
odic system and have detailed the procedure to determine 
the functions w n (r) for a general fc-point sampling of the 
BZ. Since we have in mind applications to large and dis- 
ordered systems, we restrict ourselves here to the case of 
T-point only sampling of the BZ. 

Our optimization procedure is closely related to that 
described in Appendix A of Ref. ||. We report explicitly 
the formulas to be used in a calculation with a cubic 
supercell of side L, which is the case of our simulation. 
The minimum spread criterion of Eq. (Q) is equivalent to 
the problem of maximizing the functional 



n = J2(\Xnn\ 2 +\Y n 



where X„ 



w„ 



(2) 



Similar definitions 



for Y mn and Z mn apply. Maximization of ft is per- 
formed using a steepest descent (SD) algorithm. We 
start the procedure by constructing the new matri- 
ces X^\ and Z^> via the unitary transforma- 
tions AW = exp(-A^)X^exp{A^) (and similarly for 

y« and ZW), where X$l = (wi? \e- l2 c x \w { n ) ) and 

Wn (r) = tpn(j) are the Kohn-Sham (KS) orbitals ob- 
tained after a conventional electronic structure calcula- 
tion, is an antihermitian matrix corresponding to a 
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finite step in the direction of the gradient of f2 with re- 
spect to all the possible unitary transformations given by 
exp(-A): AW = At (dQ/dA)^, where At is the conven- 
tional SD "time-step". The gradient dil/dA mn is given 
by the sum of [X nm (X* n - X* lm ) - X* mn (X mm - X nn )] 
and the equivalent terms with Y and Z substituted in 
place of X. The process is repeated for many SD it- 
erations up to convergence in the fl functional. The 
maximally-localized Wannier functions are then given by 
w n (r) = Hiexp(—A^)rp n (r). The coordinate x n of the 
n-th Wannier-function center (WFC) is computed using 
the formula 

x n = -^Im\n(w n \e- l ¥- x \w n ^ , (3) 

with similar definitions for y n and z n . Eq. (|^) has been 
shown by Restaa to be the correct definition of the ex- 
pectation value of the position operator for a system with 
periodic boundary conditions. The computational effort 
required in Eqs. (j^) and ^ is negligible, once the scalar 
products needed to construct the initial X(°\ and 
have been calculated. 



III. RESULTS AND DISCUSSION 

Amocnhous silicon has been one of the first systems 
studiedtH, with the Car-Parrinello molecular dynamics 
methodtll. We analyse here some selected configurations 
taken from a_molecular-dynamics simulation performed 
by Chiarottifcj. These configurations have been obtained 
by quenching from the liquid state a sample of 64 Si 
atoms, contained in a cubic supercell of side 10.86 A and 
periodically repeated in space. The total length of the 
simulation run was about 10 ps. Since our study is based 
on the local density approximation (LDA) to density- 
functional theory (DFT), each orbital is occupied twice 
and unpaired spin defects cannot be observed. However, 
such defects are expected to have a low density, as sug- 
gested by electron spin resonance experimentsEl 

In amorphous silicon most of the atoms are tetrahe- 
drally bonded (sp 3 hybridized); however different kinds 
of defects cap_.be -present and have been proposed in 
the literatureE3EJ~tZI: twofold coordinated atoms form- 
ing spinless, neutral defects; threefold coordinated atoms 
with neutral or charged dangling bonds; fourfold coordi- 
nated atoms characterized by stretched ("weak") bonds 
and by bond angles that are rather far from the tetra- 
hedral angles; and fivefold coordinated atoms ("floating 
bonds"). 

The usual analysis of the bonding properties is based 
on the coordination number, i.e. the number of atoms ly- 
ing inside a sphere of a chosen radius r c centered on the 
selected atom. Such a simplified structural analysis is 
sensitive to the value chosen for r c . More importantly, it 
is also completely blind to the electronic charge distribu- 
tion, which ought to be important to any description of 
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FIG. 1. Si-Si (dashed line) and Si-WFC (solid line) pair 
correlation functions. The detailed structure, in the range 
0.0-1.5 A, is shown in the inset. The data have been obtained 
by averaging over 17 configurations of the MD simulation. 

chemical bonding. Analysis of the full charge distribution 
and bonding in terms of the Wannier functions is rather 
complex. However, as we show below, the knowledge of 
the positions of the WFCs suffices to capture most of the 
chemistry of the system and to identify its defects. 

The ionic structure of a disordered system is best de- 
scribed in terms of correlation functions. We treat here 
the centers of the localized Wannier functions as a sec- 
ond species of classical particles, and we regard amor- 
phous silicon as a statistical assembly of these two kinds 
of particles, the Si ions and the WFCs. In Fig. | we 
show g(r), the standard Si-Si pair correlation function, 
together with g w (r), the Si-WFC pair correlation func- 
tion. As can be seen, g(r) and g w (r) exhibit strong peaks 
at ~ 2.4 and at ~ 1.2 A, respectively, thus indicating that 
the electronic charge is mostly localized in the middle of 
the Si-Si bonds, as expected in a covalent-bonded system. 
The spread in the peaks indicates the disorder-induced 
strain in the Si-Si bond. However, g w (r) shows some 
structure also for values of r in the range 0.5-1.0 A (see 
inset). This means that a few WFCs are anomalous in 
being very close to the Si ions. This behaviour represents 
a clear indication of the presence of defects in the system. 

If the coordination number is computed by integration 
of g{r) up to r c — 2.80 A (the position of the first min- 
imum), we find that, on average, 96.5 % of the Si ions 
are fourfold coordinated, while 3.5 % are fivefold coor- 
dinated. Therefore, according to an analysis restricted 
to the ionic coordinates alone, all defects in the sys- 
tem are identified as fivefold-coordinated Si ions. This 
is in agreement with the findings of previous ab-initio 
simulationsE3. However, when we choose a coordination 
criterion based on the Wannier function representation, 
we get rather different results. We will say that a bond 
exists between two Si ions when they share a common 
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WFC located within r w = 1.75 A of each ion, r w being 
the position of the first minimum of g w {r)- With this 
convention, we find that 97.5 % of the Si ions are four- 
fold bonded; of the remaining ions, only ~ 0.6 % have 
5 bonds, while the others are more or less equally sub- 
divided into twofold-bonded and threefold-bonded ions. 
Therefore, although the total density of defective atoms 
that we obtain is similar to that coming from the bare 
coordination analysis, the nature of the defects appears 
to be different. 

This fact is best illustrated by looking at the bonds 
formed among a few ions, in selected configurations of 
the molecular-dynamics simulation. In Fig. |^(a) ion A 
is fivefold coordinated and has 5 bonds, while ion B is 
fourfold coordinated but has 3 bonds only. In fact no 
WFC is found between ion B and ion C. Notice that the 
bond between ion A and ion B is somehow anomalous 
since the distance from the corresponding WFC to ion 
B is considerably smaller than to ion A (0.87 and 1.56 
A respectively), and the A-B bond appears to be dis- 
torted. As the ions move, the electronic configuration 
also changes, and in fact we find that after about 10 ps 
the WFC located between ion A and ion B comes still 
closer to ion B (see Fig. |](b)). The distance is reduced 
to 0.57 A, in such a way that the A-B bond is broken or 
at least severely weakened. In this configuration, accord- 
ing to our criterion, ion A is fourfold bonded, while ion 
B has only 2 bonds. Interestingly, the twofold-bonded 
atom was proposed by AdlerEJ as the lowest-energy de- 
fect in amorphous silicon. The defect we observe in Fig. 
||(b) probably represents a transient state in which ion B 
breaks the bond with ion A and tries to form a new bond 
with a different nearest-neighbour, possibly ion C. This 
conclusion is supported by the fact that the direction of 
the vector connecting ion B to the anomalous WFC is 
intermediate between the B-A and B-C directions. Fur- 
ther confirmation comes from inspection of the isosurface 
of the electronic charge distribution associated with the 
anomalous Wannier function. Notice that the density 
profile of this Wannier function is different from the one 
associated with a normal covalent bond, as shown in Fig. 
^|(a). To determine whether a new bond with ion C is re- 
ally formed or whether the configuration of ion B in Fig. 
^(b) is stable would require a longer simulation, which is 
beyond the scope of the present work. 

We stress again that the interesting transformation in 
the bonding properties of ions A and B, which we have 
described using the Wannier function analysis, cannot be 
detected using the simple coordination number criterion. 
In fact, in the configuration of Fig. ||(b), the coordina- 
tion numbers of ion A and B are the same as in the 
configuration of Fig. ||(a) , although the A-B distance is 
significantly increased from 2.38 to 2.59 A. 

The Wannier function analysis allows a description of 
the electronic charge distribution in terms of well-defined, 
localized functions; the clear representation of the bond- 
ing properties of the system based on the positions in 
real space of the WFCs can then be followed by a more 
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FIG. 2. Snapshots of 2 different configurations of the MD 
simulation of amorphous silicon. Large grey balls denote Si 
ions, while small black balls denote WFCs. A, B and C label 
the Si ions, whose bonding properties are discussed in the 
text. For clarity only the Si ions and WFCs lying within 4 
A of ion B are shown. We also plot the isosurface densities 
Pn(r) = |u> n (r)| 2 corresponding to a normal covalent Wannier 
function in (a), and to an anomalous Wannier function close 
to ion B in (b). Thin lines in (b) indicate directions of possible 
bonds of ion B with ions A and C. 

quantitative study. For instance, one can relate specific 
features of the electronic density of states to a particular 
Wannier function w n (r), defining a "projected density of 
states" 



N n (E) = I Kl^m> \ 2 S(E - E m ) 



(4) 



where V'm and E m are the KS eigenvectors and eigen- 
values. As can be seen from Fig. ||, the N n {E) function 
associated to the anomalous WFC of Fig. ^(b) exhibits a 
strong peak which matches the peak in the electronic den- 
sity of states located above the valence-band edge. This 
peak can be associated with the highest-energy occupied 
KS state. We have therefore a clear example in which 
the electronic states of the structural defects present in 
our sample are introduced into the energy gap. 

Another important advantage of the use of the Wannier 
functions is that a precise calculation of the localization 
degree of the electronic charge is possible, in contrast 
with previous approximate estimateslilllO. In fact one 
can easily compute the quantity 



(5) 



which corresponds to the spread in real space of the Wan- 
nier function w n (r). Considering again the anomalous 
WFC of Fig. |(b), we find that a n = 1.94 A, a value 
significantly larger than the average spread obtained by 
considering all the WFCs, a = 1.38 A. We find therefore 
that the Wannier function associated to the defect state 
is less localized that the Wannier functions associated 
to normal states. This interesting result confirms, in a 
quantitative way, the-qualitative observation of Fedders, 
Drabold and KlemmliZI who pointed out that the defect 
states can be far less localized than expected. 
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FIG. 3. Histogram of the electronic density of states (thin 
line, the dashed portion indicating the conduction band) com- 
pared with the projected density of states N n (E) (thick line) 
corresponding to the anomalous WFC of Fig. 2(b) (see text 
for definition) . The histogram has relatively large fluctuations 
because only a single ionic configuration is considered. The 
curve N n (E) has been smoothed with a gaussian broadening. 

IV. CONCLUSIONS 

In conclusion, we have described a novel application 
of the Wannier function analysis to a disordered system. 
We have shown that a simple geometrical analysis of the 
positions of the WFCs is already sufficient to extract use- 
ful information about the bonding properties of the ions, 
particularly in those interesting defective configurations 
which, using traditional approaches, can only be studied 
in a very crude way. In addition, since in our method 
the localized Wannier functions are explicitly available, 
a quantitative analysis, which allows us to estimate accu- 
rately the degree of localization of the electronic charge, 
is possible. It also allows us to clarify the nature of the 
relevant features in the electronic density of states. 

An additional advantage of the current approach of an- 
alyzing the dynamics of the system of ions plus WFCs is 
that it should be possible to extract information about 
the dipolar fluctuations, i.e. about the dielectric response 
function x(k, u>). In fact the present approach is closely 
related to that of Ref. |l8| where, however, only the k = 
response relevant to optical probes is computed, so that 
only the macroscopic polarization is needed at each time 
step. The current approach opens the possibility of ex- 
tracting some local (i.e. fc-dependent) information as 
well. 
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